A/ /is A C/S--/S9, //7 


NASA Contractor Report 159117 


NASA-CR-159117 

1979 00 00^20 


THE ALTERATION OF PROFILE ANALYSIS TO 
ACCOMMODATE TESTING FUNCTIONS 


Raymond H. Myers 


RAYMOND H. MYERS 
Blacksburg, Virginia 24060 


NASA Purchase Order L-72497A 
July 1979 



UNGLEY RESEARCH CEUTIR 
library, NASA 
{■lA^iSX'^LL 


rvi/\sA 

National Aeronautics and 
Space Administration 

Langley Research Center 

Hampton. Virginia 23665 




I 


THE ALTERATION OF PROFILE ANALYSIS TO 
ACCOMMODATE TESTING FUNCTIONS 


The purpose of this task is to develop a methodology for testing for 
differences between several pilot describing functions, where the data 
points represent averages at various frequencies. Typically a functions 
represents a treatment or treatment combinations which are to be compared. 
Data are taken in replicates and an output is measured for each frequency. 
It is not unusual to experience 4-8 functions representing a like number 
of experimental conditions in a real-time piloted aircraft simulation. 
Typically, as many as 8-16 frequencies may be Involved. 

Consider the following plot for the special case of two functions. 


function 

value 

(degrees) 


^22 ^^24 



frequencies , rad/sec 


Here, the represents the population mean for the j frequency under 
the i^^ function. The problem is to determine whether or not, simultaneously, 
the means for function 1 are equal to the means for function 2 , across 
frequencies. The extension to more than two functions is also of con- 
siderable interest. 


Specifically the goals of the task are as follows: 

(i) Develop and describe the methodology for testing for differences 


between functions. 
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(ii) Determine how to approach the problem of assessing the power of 
the test. Supply charts or tables for power. 

(iii) Use the power results to recommend how to design such experiments; 
for example determine an effective number of frequencies; in addition, 
what is the maximtun number of functions that will allow a reasonably 
sensitive test. Also what is a reasonable number of replicates? 

(iv) Discuss software considerations. 

Basic Assumptions for the Experiment 

In this section we shall discuss the distributional and data structure 
proposed for the test on k functions (k>2) . Consider initially two func- 
tions , in the context of Figure I. Suppose that we have two random vectors 



where x, is the scalar output for the i^^ function and the p*"^ frequency. 
iP 

These are basic output measurements. It is assvimed that each of the 
vectors follows a multivariate normal distribution with common variance- 
covariance matrix E, the latter a txt matrix. The practical implication 
here is that within each function the observations are correlated and the 
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correlation structure is the same for each of the two functions. For the 
first function n^^ independent vectors are taken and for the second function 
1 I 2 vectors are observed. 

The extension to more than two functions is obvious. It is assumed 
that there are k^ultidimensional vectors 2i2****’^ * 


*11 


*21 


*kl' 

*12 


*22 


*k2 

• 

« 

• 


• 

• 

• 

... 

• 

. 

• 

*lt 


*2t 


*kt 


The assumption is made that is multivariate normal with mean 
(i=l,2, . . . ,k) and variance-covariance matrix Z. 


Hypothesis 


Consider again the case of two functions illustrated in Figure 
A reasonable hypothesis that accomplishes the goals outlined here, is 
given by 




= ^21 
^^12 " ^22 
^13 “ ^23 
^14 " ^24 


( 1 ) 


or, in terms of vectors. 


*0- % ■ JI2 ’ 
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where v. is the mean vector for the i function. One can easily see 

— X 

that this is the hypothesis that the experimenter needs to test. Namely, 
the null hypothesis states that the mean value of the two functions are 
equivalent ^ each frequency . In the general case of k functions and t 
frequencies the hypothesis is 

"o' 'll =H 2 '^3 ° ••• 

where each vector is t-dimensional. We need a test procedure for testing 
the joint statement given in (2)* In addition, we need guidelines on 
the experimental design that will adequately give evidence in favor of the 
alternate hypothesis when it is appropriate. 

Profile Analysis 

The problem posed here appears to be very similar to Profile Analysis 
[2] which is a multivariate technique in which a "repeated measures" 
design is used and the treatments are broken into groups. Basically, 
measurements on individuals (or individual items) are conducted for each 
of say t treatments, with the observations correlated from treatment to 
treatment. The purpose of the experiment is to determine if there are 
differences in treatment means. When the treatments are put into an 
additional classification variable called "groups" then the procedure 
becomes what is called Profile Analysis . Consider, again Figure I and 
allow the functions to take the role of the groups and frequencies take 
the role of the treatments. This would appear to be an application of 


Profile Analysis. 
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In a Profile Analysis the test on treatments is not valid unless no 
interaction exists between treatments and groups* Thus* initially a test 
for parallelism (no interaction) is usually conducted before the test on 
treatments is attempted. The hypothesis of the test on parallelism is 
given by 

^11 " *^12 ” ^21 " ^22 
^13 ” ^12 “ ^23 ” ^22 
^lA ~ ^13 “ ^24 “ ^^23 

The test involves a Hotelling’s T^ [2] and is quite easy to extend to 
more than two groups. 

Alteration of Profile Analysis to Accommodate 
Testing Functions 

One can easily see (easier to visualize in the case of two functions 
or groups) that if groups and interactions between groups and frequencies 
(functions and frequencies) are not statistically significant, one can 
interpret this as implying that functions are not signxficantly different 
point by point or frequency by frequency . One mode of verfication of this 
is to note that if two functions are parallel between frequencies the only 
way that the two function averages can be the same is for the two functions 
to coincide. On the other hand, function averages might be very close 
together but because of non-parallelism (interaction) the two functions 
may be far from coincidental. Thus neither a significance test on groups 
(function) averages or interaction will suffice for our purposes. Rather 
it is the joint hypothesis of the two that is relevant. 
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Searching into the structure of the problem more closely suggests 
that the condition of coincidental functions we have described here is 
identical to the hypothesis in (1) or (2) for the general case. Thus it 
can be said that we have a profile analysis in which we are attempting 
to simultaneously detect parallelism and group (function) effects but we 
are doing it with a single hypothesis and thus a single test and not with 
two tests. 


Test Statistic and Procedure 

For the hypothesis in (1) , we are simply testing the equality of 
two mean vectors. For the j function, n^ Independent vectors (replicates) 
are obtained and we have two averages which are joint estimates of 

the mean vectors and y ^2 (t-dimensional) . The sample data is used to 
compute an estimate of the variance-covariance E. Thus we have for the 
first function, n^^ independent vectors, and ri 2 independent vectors for 
the second function. The estimate of Z is found by obtaining sample 
variances and covariances for each function and pooling over functions. 

Thus the estimates have n^^ + TI 2 ~ 2 degrees of freedom. The test statistic 
is given by [2] 

T^ = (x^-X2)S ^(X]^“X2^ 


S = £(— + — ) 


n. 




where 
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The statistic in equation (3) follows Hotelling’s distribution with 

n + n- - 2 degrees of freedom. If the same number of replicates are 
X 2 

taken for each function, say n, then of course n^^ + U 2 ~ 2 is replaced 
by 2(n-l). The statistic 


(nj+n2-2)t 


(4) 


follows a F distribution with t and n^^ + n 2 - t - 1 degrees of freedom. 
This fact allows us to handle the power quite easily. 

The extension to more than two functions is quite easily done. 
Consider the hypothesis of equation (2) , which, written in expanded form 
is given by 

“o= >'ll " ’'21 ° “31 ' ••• '“kl 
“l2 ” “22 “ “32 “ ■ ■ ■ “ “k2 


“it ■ “2t ■ “3t 


= y 


kt 


Once again , there are vectors of averages , X 2 » • • • » 2^ » t-dimensional . 
This data is used to compute a pooled estimate of the variance-covariance 
matrix I. In addition a txt matrix (hypothesis matrix) is given by 


1 1 

h = E - T . T . - . 

rs nj rj js N r s 


Here T. is the sum of the data on the r^^ frequency for function j. 


th 


N = E n. . G is the grand total of all observations on the r 

j=l ^ 

frequency. The test statistic is the largest eigenvalue of the matrix. 
HE“^ , where E is the error sum of squares matrix , which is given by 
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E = E • d.f. 

where d.f. are the degrees of freedom for estimating the variance- 

k 

covariance matrix (d.f. = Z n . - k) . 

i=l 

Examples of the computation for both the special case and the general 
case will be given in a later section. Charts for the use of the largest 
root statistic can be found in Morrison [2]. A computer program for the 
computation of the test statistic will be found in the Appendix. 


Power of Tests 


The power associated with the case of two functions is quite easy to 
handle from a theoretical point of view. In a practical sense there are 
some difficulties that are not insurmountable. The test statistic in 
(3) follows F 4 _ _f 1 under the hypothesis of equality of means 
stipulated in H^. One can immediately notice that the denominator 
degrees of freedom automatically imply an experimental restriction, 
namely that n^ + n 2 > (t+1) . The test procedure, of course, involves 
rejection of if the calculated test statistic in (4) exceeds the 
upper tail percentage point of the F-dlstribution. If the test 

statistic follows the non central F-distributlon [1] with degrees of 
freedom t and n^ + n 2 - t - 1, and with non centrality parameter 


n n 
= 12 
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Thus for a specific difference between the mean vectors and a 

specific variance covariance matrix , the power is given by the prob- 
ability 

' 'OL’ctM and is 

the upper percentage point of the F distribution« Thus the probability 
distribution of the non central F is needed. This distribution is 
discussed in detail in Graybill [1]. Power charts are offered in 
Morrison [2] and are quite easy to use. An explanation of their use is 
given in the text. 

In order to compute the power one must know Z, or at least have an 
estimate of it. For the study described here, hypothetical values for 
the variance-covariance matrices are used, these values generally taken 
from previous studies at NASA Langley. A power study is described in 
the following section. This study was done to arrive at appropriate 
values of sample sizes (number of replicates) and number of frequencies 
that provide good power for specific values (in terms of order of magnitude) 
of the differences (_Uj^-p^2) being detected between the functions. 

The power computation for the test for more than two functions is 
extremely difficult and is wrought with practical difficulties that 
actually prohibit its use. Approximations are available [3]. However, 
the approximations involve many queintities that must be estimated from 
the data. The non central F probabilities are monotonic increasing 
functions of the non centrality parameter which, if one investigates 
the structure of 6^, is very pleasing. That is, 6^ represents "how 
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false" the null hypothesis is before we would wish to reject with a 
high probability, and as grows large we would expect the power to 
increase monotonically. However, the structure of the approximations 
of the power of the "largest root test" is such that several non centrality 
parameters appear and the power is not a monotonic function of them which 
suggests that for at least some cases the approximation is not particularly 

good. 

It would seem reasonable then that sample sizes and other recommenda™ 
tions should be done for testing two functions at a time . In addition, 
the user can assess his own power in a particular case by computing the 
power through the use of tables of the non central F . This will be 
reasonable as long as one does not study too many functions simultaneously . 
This does not mean that one must restrict the number of functions in the 
test, but rather in a power assessment it would only be accurate if the 
number of functions were kept low , say four or less . It should be 
emphazised here that the value, t, the number of frequencies, n^, the 
number of replicates are more crucial than k, the number of functions. 


Example of Test Computation 

In this section we present some hypothetical data and numerical 
examples of the computation for cases of two functions and three functions. 
The output here is not what NASA Langley personnel has at its disposal 
but it does represent a representative type format for the output. For 
the case of two functions there are t = 5 frequencies and n^^ = n 2 = 
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First, all univariate F-statistics are given comparing function values 

A. 

for separate frequencies. In addition, the 21 and H matrices are given. 

The Hotelling’s value is given for the case of two functions and the 
largest eigenvalue of E”^H is given for the case of three functions. 

In both cases, the test statistic is significant. These examples are 
shown in the Appendix. 

Recommendation of Experimental Design 
Through Power Study 

An extensive Power Study was made in which the power of the test was 
computed for varying number of frequencies , number of replicates , correla- 
tion structure, and assigned differences between the function means. A 
set of plots are provided in order to give the reader an indication of 
what magnitude of power can be expected for various combinations of the 
parameters and to devise recommendations of what are reasonable values 
for the number of replicates and frequencies, given a particular correla- 
tion structure. It should be emphasized that regardless of power con- 
siderations there is a certain restriction that must hold. For the case 
of two functions, if we call n the number of replicates for each frequency- 
function combination, the denominator degrees of freedom for the F 
statistic will be non-positive unless 

2n > t - 1 (5) 

Thus equation (5) represents a necessary restriction in the case of two 
functions. In the general case it would be adviseable to keep the same 
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restriction, i.e., that given by equation (5) because it is unlikely that 
the quality of the test (i.e., power) will be good at all unless n 
satisfies equation (5). In addition since is is recommended to compute 
the power on the basis of Hotelling's (tests on two functions) one 
cannot assess the power unless equation (5) holds. 

We attempt here to give some rationale regarding our power study. 

The most difficult variable to cope with was the covariance structure. 

Keep in mind that the correlations represent association between observa- 
tions from one frequency to another within functions. The following 
covariance or correlation structures were studied: 

(a) Low correlation, correlation constant (0.1, 0.2) 

(b) High correlation, (0.6-0. 9) 

(c) Mixed correlation between high and low (0.1-0.95) . 

In addition, the power (through the noncentrality parameter) is a function 
of ~ ^ 2 * difference between the means of the two functions. Here, 
we are assuming that there is a typical difference which is a proportion 
of^ £, where a is the same at each frequency. Proportions of 0.3, 0.5, 

0.7, and 1.0 are used in this study. We feel that it is reasonable to 
assume that 0.7a - 1.0a is a reasonable range in the difference between 
the two functions which one would like to detect. It is clear from the 
plots that for any value less than 0.7o, experimentation is much too 
demanding if one wishes to obtain high power. The significance level of 
the test was fixed at a = 0.10. Plots of power against frequency are 
shown in the Appendix for various correlation structures. Curves are 
shown for various proportions of a which one is attempting to detect. 
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The purpose of the plots Is to attempt to show the role of t, the number 
of frequencies and n, the number of replicates, on the power. The plots 
indicate the following: 

(a) All correlations equal - low correlations (0.1, 0.2). The 
role of t, the number of frequencies is displayed dramatically in the 
plots shown. As t grows large, the power is reduced, and for t > 10 the 
only condition under which power is moderately good is when n > 10. 

Thus one can assume that the test will be ineffective unless t < 10; 
but if t must exceed 10, the number of replicates should be 10-12. 

(b) High correlation (0.7-0. 9). Here it is very clear that despite 
the number of frequencies, in order to obtain any effectiveness in the 
test, at least 12 replicates should be used. Of course, it is probably 
an unusual experimental situation for which all correlations would be 
high. 

(c) Correlations mixed between high and low. This condition likely 
represents the most typical experimental condition that one confronts. 
Again, however, it appears as if small numbers of replicates will not 
give power values that are acceptable if one desires a quality test. 

For 8-10 frequencies, n = 10 is barely acceptable whereas for t > 10, 
a sample size of at least 12 is necessary. 

Conclusions 

A test procedure is given for testing equality of k different pilot 
describing functions. The method for handling the power for the case of 
two functions is described. It is felt that this is not an unreasonable 
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assessment of power as long as the number of functions is four or less. 

For rather typical types of correlation structures a power study is made 
in order to determine reasonable number of replicates and frequencies. 

Two other recommendations come to light, in addition to those involving 
the experimental design. These are given in the following two paragraphs. 

In using the techniques described here the NASA Langley personnel 
should be very careful to sample to determine what correlation category 
prevails. The power is very much dependent on the correlation. In fact, 
if possible it might be adviseable to take data Initially to determine 
the approximate correlation structure, then supplement with the required 
data according to the experimental design recommendations. 

The user should not be attached^ to a particular significance level. 
We used a = 0.10 for the power study. However, because of the inherent 
lack of power in tests such as these one should be prepared to consider 
that a test which is significant at the 0.15 or even 0.20 level is 
evidence in favor of differences between functions. This is particularly 
true if one is forced to use a large number of frequencies. 
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APPENDIX 


Computer Program and Writeup. 

There are three subroutines supplied in the appendix for calculating 
various multivariate statistics for many hypothesis testing situations. 

Some notation: 

H — The hypothesis sum of squares matrix* 

E — The error sum of squares matrix* 

NP - The dimension of H and E* 

NUH - The hypothesis degrees of freedom. 

NUE - The error degrees of freedom. 

Xhe first subroutine ROY requires as input: 

H, E - Square matrices of order NP stored columnwise in the vectors 

H and E respectively. 

NUH, NUE - Degrees of freedom are also supplied by the calling 
program. 

RTS - A work vector of length NP. 

N - A work vector of length NP*NP . 

Output: 

RTS - Contains the ordered eigenvalues of HE ^ . 

theta = (-^v§2L_) „here Xmax is the largest eigenvalue of HE"^. 
'1+Xmax 

LS, RM, RN are S, M, N respectively for use in Heck's charts. 

The second subroutine NILKS needs as input H, E, NP, NUH, and NUE as 
previously described, where H and E are square NP*NP matrices stored 
columnwise in a vector of the same name. 



A-2 


Output : 

IfI 

ui = — J — ■ — for comparison with Wall's Tables. 

E+H| 

W2 = -(NUE-(NP-NU1H-1)/2)*LN(W1) for comparison with chi-square 
tables, in large samples where W2 has NP*NUH degrees of 
freedom. 

W3 = F-ratio having INUM numerator and IDENOM denominator degrees 
of freedom when and exact F-test is available otherwise 
W3=0 and INUM=IDENOM=0 . 

The third subroutine HOTEL can calculate M T^-statistics . The 
necessary input includes 

X - a vector of length M*NP containing M sets of mean vectors 
of length NP stored columnwise. 

S - Sigma - where X(J) , (J=1,M) has covariance-matrix 
(l/A(J))*Sigma. 

A(J) - A vector of length n containing . (where the n are 
the number of replicates in the i^h function) , stored 
columnwise. 

T2 - is a work vector of length M. 

Output: 

S: S -inverse 

T2 - contains the M values of Hotelling's T^ for direct comparison 
with tables 

lEH - comes from matrix inversion routine 

= 0 No error occurred during inversion 

= -1 No result because of wrong input parameter n or the 
matrix is not positive definite 
= k if there is a loss of significance. 
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SUHPCHJUnE «0Y(H,F. r MUH,NUfc#WTS^ THFTA,LS.RMf RN) 

SUHRlJUT t Ut S wtrJUlREO; ONWUOT A.N»D OElGf.H SUPPLIED 4Y USER. 

'MATRICES h DUE TU nyPoTHtSlS ANO t OlJE TO fcPROR, Olr^ENSION NP, AND 

DEGREES (»F FREfcOO'*'^, nUh FUR HYPCiTHtSIS AM) UiJt FUR ERROR/ ARE SUPPLIED IN 
T*-F, call TUG PWOURiv^ 

H A\0 F. APF SUUARt f^ATRlCES OF C»R0ER NP STORED COLUN NU I SF . 

RTS A.f;|) auRH vectors of DIMENSIONS NP AND NP^NP/ PE SPE. C T I VF L Y . 

VAI ijtS PEIUPNED AkF The vector RlS OF riROKRFO LATFNT ROOTS UF H^E> I N VERSE / 

Tir. T Ar (LAkGE ST K •. Ji I T ) / ( 1 , + 1. A PGF S T ROO T ) / A!vO * P AS AMF 1 E PS LS# RM, AE.'D RM FOR 

RFFERK.JLE To HhC^'S Charts (DROP THE. FIRST LETTER OE EACH). 

DIVEUS1«)N H( 1 ) ,E( n #rt( n /RISC 1 ) 

DOUBLE PREClSIflU H » E , R I S r T HE T A / X , a 
I - IS, IS/PO 

IS LS = \Ur. 

GO Ti: /:?S 
?0 .S = '^P 

?S R'*=CFL OAT ( 1 AhS(NOH-UP)-1 ) )/^?. 

J'Jr(FLUAl (Nul-NP-1 ) )/P, 
call ON^urjl (T.PrH.Ef RTS,‘^i) 

IHETA=RTS(1)/(1.UU*PISC1)) 

RETURT4 

END 


St'HROUT INF iVlLKS(H, E/NP/NUH,NUE/in1 / ^j2/ w3 / 1 NIJM , IDENO^'W IEw) ^ 

SJRROUTINFS REUUIRED: DMFSD FROM SSP LIPRAry. w 

'lATwlCES M OLIF TO HYPOTHESIS AND F DUE TO FPPOP, DIMENSION NP, AND 
DEGREES OF EKEEDUM, NUH FOR HYPOTHESIS AND NUE FOR ERROR/ ARF, PROVIDED 
BY THE calling program, 

H AND E AkE stored C OUJMN-rt 1 SE . A S UPPER TRIANGULAR MATRICES.' 

OUTPUT INCLUDES: 

-I =UE I (F )/OE T (E fH) FOR COMPARISON WITH vJALL'S TAbLESr 
.•4P = -(riUE-(;:P-NUH+n/^)*LN(rM ) FUP Ci^^PAPISiJH wTTH CMl-SQUARF HAVING 
NP*\Uh degrees of freedom in large SA-'^plKS- CLN(.) is NATURAL LOG). 

/i3 = F-RATICJ HAYli.G luUM NJME.RA.TOR ATjO 1?}F nO‘^ 0EN0'‘*' I N A T UR DEGREES 
OF FI^FEOOM /jhE.M exact F-TESTS AkE AVAIlaHLE: OTHER^^ISE /«3 = 6.'0/ 

Ii>jU‘'» = 0/ A>jO lOENQMsO. 

dimension M(1 )/t(l) 

OOiJHLE PRECISION h , E # T / 1 r ft 2 / w 3 / , |)E TE / OF T T / PROD 

TOL=. 000001 
_ = \p»(MP*n/2 
Oj 20 I=l/L 
T=h( 1) 4h( I ) 

20 H( n = T 

call OVFSTXt /NPr TOL/ lEE) 
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CALL O^FSO(h,NPrTOLf lEH) 

= l *1)0 

LI =1 
LP = ^ 

)n 1^3 ,NP 

(L n *pwuu 

Ll=Ll*L^> 

2S l;> = l r'* 1 

D t I r = ►" H U 0 * * d 
•^^U(J=l ,D<» 

Ll =1 
LP = P 

;j 50 I = 

->^ru! = »* ( L 1 J 
-1 =l 1 ♦L? 
iO L^=L^H 

.■): T T=PkOD**? 

M rOK re /iJb TT 

/.virOhLUAT ( -JUF )- I OF LO A T ( NP - SiUH + 1 ) )/?.00 
aV z - v\ * 0 L U 0 ( w 1 ) 

IF (MJH.E.J. n GO TO SO 
I*" ('UJH.hO.^) GO TO hO 
I - (’-►'.K J. 1 ) GO TO 70 
irc\^,e.o.^) GO 10 (U) 

S3 T(1 ‘?0 

so I = 

I 3h vi(j -^zoue +\UH-\P 

M = ( I .00-A n *DFLOA r ( IDENOM) j ♦okLOAT ( INOW!) ) 

j j TO 150 

hO 1 .'JJ‘-*z;>*NP 

I OF M(|V = P* (,vjijt + MUH-tMP-l ) 

= | 1 ) ) *OFLOAT ( IDENOM ) / ( OSQW T ( w 1 ) ★DFLO A T C I OUM ) ) 

7 0 i N ' J ■•T Z \ i. I H 
1 ?F ‘jOv.z'nIUF 

:'/5 = (1 , 00-iv n *OFLOAT (IDENOV) / i ♦OFLOAT ( ) 

.•3 r*i 15 0 
80 I OWzP 

I (i^ot-l ) 

1 .oo-nstjlvr (|«1 ) ) *OFLOAT ( IDtNOM)/(OSfJKl (ml) *))FLOA) ( INU*^) ) 

.♦I I TO 150 
A 5 = 0.t)t) 

I 'JJV = 0 
I f)LM‘^’ = 0 

150 IF (ItF GO TO I 5S 

G3 TO US 

15S IF(IFh,Mj.0J GO TO UO 
GO TO luS 
un ]E/. = o 

G) TO ISO 

los IEazI 
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DATE 03/51/7A 


PAUb* a 


IbO -^ETUPN 
FNO 

S^ihkCjUT I%F. HOTEL (Sr Ar *>JP» X , A , T<>r U.H) 

SJt^wiMIT 1 NLS WfcOUl^^bl): OSnJV IN S3P LI'^MAPY 

S IS A SV'^^tTWlC ^'ATMlx OF HPOtW NP STOPtU COLll'^N -W 1 S£ IN UPPER 
T^VI A\{, iLAh bfIKM* 

A A- i) W AkL vhC. »H^^S nf length 
I; s a . u^-^ v't-cnw np LENGTH .JP, 

ASSir’b VfCTuwS X(J) To HAVE COVARIANCE MATRICES M . / A ( J ) ) *S I GM A , FOR 
J = lr^^ TmL vaTRIX S Ai'JO VfcCTUPS XU) ARC SUPPLIED BY Tht CALLING PROGRAM 
AS ^*.ELL AS The OIMFmSTON NP AND M. THE VECTORS X(J) AWE STORED COLUMN- 
^ilSE IN the VbClOP X OF LENGTH NP*M, 

0\ »-FTijW'i S IS PtPLACF.D Hy SINVFRSF, EACH x(J) IS REPLACED BY THE 
CriRRt SPO*jDIN(; coefficients fJF A LI'viEAR 0 t SCR 1M I N ANT FUNCTlOu, ANO 
TmE VFCTur TP CftNTAIMS VALUES OF THE HOTELLING'S T-SQUARE STATISTICS 
F.TR I)IRF:CT CU-IPARISUN rtJTH SPECIAL AID TABLES (SFF JENSEN AND HlTwE r 1967) 


?0 

PS 


30 


0 I MEN'S UiN S ( I ) r X ( n r A ( n , T P ( 1 ) r ( 1 ) 
DOUHLE PRECISION S , X , A , T P r « r StJM , SUM 1 
iaL = .()onooi 

CALI t'SlNVCSrNPr TOLr lEH) 

LP=l 


3 = wP 


DJ SO 1 = WM 


S 1 .*Z0 . :)(» 


LI = 1 

DO PS J = l.PrL3 
DO PS K=LPrJ 
5J*^l=<CK)*S(Ll)»x(J) 
TF ( J-K ) ISrPO, IS 

S / •' 1 1 ■ 1 c» , i)<) 

SJ' =SU-’ + SU-^l 
L1=LM1 

TP(I TzSlJvUAd) 

LP = LP>'*P 
L3=L3*NP 
. = 0 


Lr = ur 

DO SO I=lrM 

L^=L*NP 

OJ an J=LlrLP 

SJHzO.DO 

00 3S K = LI rLP 

LJ=J-L3 

L-'sK-Li 

call LOC(LJrLK.J*^rNP,NPrl) 
SUM=SUV+S(JK)^X(K) 



o-iTf 


3S a(LJ)=SUv 
UO COM PiUt 
JJU as JJ = L1 
JU = JJ~L 5 
as < r j J } ( iL ) 

L = L ♦ I 
LI =L1 
SO 

^ u*' s 

bND 

SUS^OUT Pit DfviKnnT (M, XL r X ) 

D! ^^t'.Sins A( 1 ) ,M( 1 ) , XL ( 1) f X n ) 
OClUHLb PNetibllUJ ArH,XL»X,SUMV 

< = I 

DJ 100 J=^,M 
. = '•* (J-1 ) 

00 100 1 = WJ 
u = L*l 

< r ♦ 1 

lOO )=rUL) 

•' ‘ z < » 

CALL Ot P;tN(hf X, 

L = 0 

DO no j=nM 

L = L + J 

IP* ' . I ' ; = 1 . ■'/ ^ (^A-S( .(D)) 

d 5 ns j = nw 
DO ns i = n^ 

< = K ♦ 1 

ns r<(-0=>.(K)*XL( J) 

DO 1^0 l=l,Ni 
\?z:0 

DO 1?0 J=1,M 

\'i=‘^*M i-n 

L = '"^{ J-1 )^*I 

(L)=w.00 
DO 1<^0 K = 1,M 
VI =‘;l ♦ 1 
V?=N?+1 

1^0 X ( L ) =* U ) ♦^KNl) A A (N;^ ) 

L-“ 

O') MO J = 1#M 
DO MO i = nj 

V 1 = 1 -A' 

■ ( J-i ) 

L = L * 5 
A(L)=''-P‘0 
DO 130 Kzn.'^ 

VI =1>I1 f M 


PAUL S 


NPno 3/0 
Nwno 3H0 
Npiia aso 
MPOh hlO 
Npon spu 
rjkoo ^^30 
f)ij(ui h/i n 
riPDfi sso 

NWIU.I hh,i) 

NPOO h70 
NKOfJ 710 
NWUn 7?0 
NPOO 770 
NPOn 780 
ivpun 7^0 

MPOO 80U > 

NPOO 810 
NPUO 8^0 
NRiJU 8 50 

N S n G 8 a 0 
NPOU 8Sn 
fNjWOO 8^0 
Nwf)n ^00 
NPOl) 910 
NkOO 980 
NWGG 930 
rj R 0 n 9 a 0 

N9GG 9S0 
N'kUll 98 0 
NWOO 970 
NkGG 980 
NPlll) 990 
Nk(j01on0 
vkuniolo 

NROUl 080 
(VRUil 10 30 

, . k ( . ( n n -4 0 
ukCifil OSO 
NROU1O80 

NROOIO/O j‘ 



I.* 


DATE 03/^1/78 


PAGE 6 


no A(L)=ACL) 

:aLL *j-irJb''itA,x,vi,vv) 
L = t‘ 

OJ luo I = l,v, 

L = L^I 

luo XL( I )=A(L) 

!M’ iSn I = WM 
NiPsO 

DJ no jsi ,vi 
\) = I-^• 

L= '*(J-1 )♦! 

i‘(L)=C'.iJO 

■ : n ■ = ! . 

\i 1 r ‘ 1 ♦ 

= 1 

ISO a(L)=A(L)*M(M)#X(Nf>) 
LfO 

" ' 1 •'> J ■ 1 . ' 

S.KV = n,no 

o:^ 170 i = 

L = L-H 

170 S.;^*V = SU'^VfA(LnA(L) 
17S S.»v^v = i'S'J^7 (Su.'-W) 

OJ tMh I=J ,v 

IPO X (K ) =A (K l/SU- V 
^ n • j • ' • 


NROD1080 
NPOOl ORO 
NPOni J 30 
NKOin 1 UO 
iNPOUl 150 
NPoul no 
NROOl 1 70 
NPoumo 

NPfUM t?/^0 
•'JPOOl ?30 
NPOf) 1 etio 
riPf)n 1 ^50 

r*piif» 1 c->/so 

f. M.. • n 70 
nKO 
^Pll^l^:>oo 
h p n n 1 M) 0 
hkUdl 310 
uPuU 1 

i I n * 0 
IMPIJO) 3U0 
dP(U' 1 3S0 
fMHfrjfil no 
»MK001370 
ivpoon^o 

dKfiO 1 300 
N ru) 0 1 u 0 0 
dwnn 1 u 1 0 
f'*PMO 1 U PO 
P j u iO 


L-y 
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• EmbFP .'jAMfc EK;F.fi 

PIGE 

eic;e 

sU'j^^nijTiME figfnj lint 

FIGE 

^^MWP(»SE EIGF 

hIGE'MVALUtS AND FIGEMVFCTUPS OF A AL SYMMETRIC tlGF 

M A T 1 X F 3 [; F 

FCGE 

US Af,f j 

call f.1c;fn( a,w,r,.^v) tJGt 

F I r F 

D t S C ►< I R T I * I i'.' u F R A k A V t I E S E 1 G F 

A • n»^lGliNAL MATRIX ( symmetric ) » DESTROYFO IN COMPUIATION. EIGF 
RtSULrA'JT eiGKu ^/ aloes are DFVtLnPED IN DIAGONAL OF ETGE 

•-MATRIX A IN UFSCFUOI'IG OR{lF R . FTGF 

R - RtSOLlA.NT MA]Ri< fiF F UJF -^\/KCTORs (STOPF!) C( 5 HjMM%ISF, PiGh 

IN SAVF vSFU!iFNCF A.S E 1 GF N V ALMF S ) EIGP 

N-ORiTFRUFmaIRICESAAUDR FI(;F 

A'v- INPUT coot §ir,F 

0 CU‘-‘RIPTF FIGE-NYALUFS AND E I GF N V T OR S FIGF 

1 COMPUTE EIGtWALiirS ONLY (R'mFFD NOT BE FlGF 

r)l.''E<vSiONhl' H(JT MUST still AP^F.AR IN CALLING FIGF 

StUUFNCF) FIGE 

REMARKS EIGF 

ORIGINAL matrix a mjst HE REAL SYM'.^P7PIC fSTORAGE Mf,DF = l) F3GL 

••'ttlRlX A cannot Hh IN THE SAME' LOCATION AS MATRIX R ETGF 

SOrtROUT If^ES and FUNCTION SUBPRIIC.R AMS RFOtiiRFO FIGE 

NONE EIGF 

E T H u i.) F I G E 

01 AGUNALI/AT lOM METHOD ORIGINATED By JACOHT AND ADAPTFO FIGF 

f^Y VON NFUMANu FOR LARGE COMPUTERS AS FOUND IN ' M A 7 HF M A T 1 C A L F 1 G E 
^FThUOS for digital COMPl.tTFPS'r F.DITFD BY A. RALSTON AND FIGE 

H.S. I'^LF, J()H\ wIlFY AfU) SONS, YORK, IROP, f.HAPTFR 7 EIGF 

ek-f: 

,E iGt 

SJHPOUTINE EIGEN(A,r,N, MV) llnp 

DlMFNSltlN A(l),R(l) EIGF 

f:h;e 

• . . .EIGE 

VERSION OF TFPl^^wnuff^MS DF.SIRTD, THE EIGE 
I r.. WEVOVEO FROM Thf UOUBLF PRECISION EIGF 

STATEMEM aHICh FuLLOaS. EIGE 

DOUBLE precision A,R, ANORM , ANR--X , T mR , X , Y , S I NX , S T NX ? » C OS X , MGE 

1 cosxp, SI Ncs , range t;iGE 

ImL 9 r*-- FROM DOrjRLE PRECISION STATEMENTS HgE 

APt*EAR]NG IN OTHER ROUTINES USED IN CUN.IUNCTTON aITH ThJS EIGF 

R u 0 T 1 E . E I G t 

The DOUBLE PRECISION VERSION uF THIS SlirtROUTTNl MUST ALSO MGF 

contain DOlJHLt PRECISION FORTRAN FUNCTIONS. SORT IN S I A T EMEN T S E I G E 


10 

ao 

30 

BO 

hO 

70 

BO 

90 

100 

no 

IPO 
I 50 
1^0 
ISO 
I HO 
170 
IBO 
190 
poo 
PI 0 
PPO 
P30 

poo 
PSO 
PftO 
P7n 
PRO 
PRO 
300 
310 
3P0 
330 
3/40 
3S0 
3B0 
370 
3A0 
390 
aoo 
ai 0 

apo 

^50 

iiUO 

tihO 

a70 

aHO 

'4Q0 

BOO 

BIO 

SPO 

B30 

S/40 

BSO 

BBO 

S70 

BBO 
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oorsnoono 
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NAMF t. lOtN 
a 0 » t 7 b 
'-!ijsr Hfc 


Ah') 7A MiJST CHA\G£() TO OSQPT 
CHANGbO TO OAHS. THt CONSTA^aT IN 
7 U 1 . 00 - 12 . 


GF.Nt-^ATL lOtNTITY '^AT^nx 


AHS IN STATfcMENT EI^E S^O 
STArtMPNT 5 SHOULD EIGE HOO 
EIOE MO 
tIGF 6?0 

Flop. 630 

EIGE 6U0 
EIOE 6bn 
F. IGE 660 


C 

C 

C 


r 

c 

c 


c 

c 

c 


S y A\r.F = 1 . OF -6 

IF ( -iv-1 ) lO.Pbf 10 
10 liw = -'\ 

(;h 20 J = 1 r*>< 

1 

00 20 i= 1 ,N 

IJ = 1<>-*1 

( IJ ) = 0 . 0 

IK(I-J) 20,lbr?0 
1 S P ( IJ 1 = I . 0 
20 COMIMUt 

CO“^PUTt initial A ^^0 final NOri^S (ANOb'‘"’ AND ANO^^X) 


2^ AufJPv.ro. 0 
DO 5b l = WN 
D ij 5 b J s I f N 
IF(I-J) 50#3b.30 

30 t A = 1 4 ( J A j ) /2 

AOi'P* SMD(iHV>fA(IA)*A(IA) 

3S Cflf T I j'lF 

IF (A‘ -JNV) IbS, 16b, UO 
iiO A 'jO'- •■ = 1 , a 1 u f ( A'.;v ) 

J( = * p /F Lii A T (n) 


1MI1ALI7F. INDICATORS AND COMPUTE THHESHOLO, THP 


INO = 0 
1 HiirANfiRV 

iiS THK = TMR/FLOAT(i^) 

SO L = 1 
bS v = 

COvpilTfc SIN A.>iO COS 


MQr (\,*M-V )/2 
LU=(L*L-L)/2 

LM = L* Vrj 

h? 1 F ( AHS ( A (LM ) ) -THR; I 30 ,6b# 6S 

6b INI) = 1 
LL=L*LU 

wvr M + 

X=n,b*(A(LL)-A(MM)) 

6M Y=-A(LM/ SQRT(A(LV)*A(LM)fX^X) 

IF ( X) 70, 7b# 7b 

7S sIbX = Y/ SORT r 2 .0* ( 1 .0+ ( S;JP T ( t • 0-Y * Y ) ) ) ) 
STNx 2 rSlNX*SINx 
7rt CUS*= bORl ( I .0-SlNX?) 

COSx 2 sCOSX*CUSX 


tIGF 670 
FIOF 660 
F. inF 6P0 
FIDE 700 
FIOF. 710 
FIGF 720 
EIOE 750 
tlGF 7ao 
ElGt 7b0 
FIGE 760 
FIGE 770 
FIGF 760 
FIGF 790 
EIGF 600 
FIGF «10 
FIGF. H20 
fcIGt 630 
EIGF Hao 
FIGF 8b0 
FIGE 660 
EIGF. 670 
FIGE 660 
FIGF. 690 
FIGF 900 
fcIGE 9iO 
FIGF P?0 
EIGF 930 
FIGE 0^10 
FTGE 9bO 
EIGF 960 
EIGE 970 
EIGF 960 
EIGF 990 
EIGE 1000 
El GF: 1010 
EIGE1020 
EIGF 1030 
ETGFl Oao 
FlGElObO 
E1GE1060 
EIGE 1070 
EIGE 1060 
EIGF 1090 
EIGE I 100 
EIGF M 1 0 
ElGFllPO 
EIGE 1 1 30 
EK;E 1100 
EIGEl IbO 
ElGtl 160 


> 

I 

so 



oonoo 


SIMCS sSI'JX*COSX 

c 

c t^OT/iTt L AMD »' 


COLU'^^-MS 


(L-1 ) 

1 v.jru M '1- n 


• 1 v.jru M '1- n 

i')il I = 1 f • J 

1^J=(] *1-1 )/2 
T ^ ( 1 -L) ^0, 11 S, HO 
HO 1 M 1 -VI) 1 1 S,90 

8S I^=l+^u 
i;u Tn 

go I ♦ I 'iJ 

I K ( 1 -L) 1 00 , 1 OSf 1 ns 

1 OO 1 1 = 1 + Ltj 

un u» 110 

lOS llzl^lu 

MO x=A(Il)*COSx-a(1M)*S1NX 

A ( iv' ) = A ( I L ) * S I N X ♦ A (I M ) A c OS X 
A( lL) = x 

ns ) 1^0M^S#1?0 

l?0 = 

1 “k::! -s4 I 

x = i<( lL^")*CuSx-P( M’9 )aS1NX 

g n ..-K ) = -v ( IL^ ) *S1 v^x + ^^ ( 1 -Ik ) *cusx 

k ( 1 1. k ) = X 

IPS CO-Nilirvi't 

X = P.OaA(LN.) aSI'SCS 

Y = A(l.L) *Cu5X?4 A(M'^)*SINXP-X 

xra (LU *Sl.^xP4 A ( '^vi) ★CUSxP + x 

A(L 1) = ( A(LL)-A(M-1) ) aSI jCSf A(LM)*(COSX?-SImx?1 
A(l L)=Y 
A( v .-.)=X 

I f- S I S HJk Cuv'PL^ T ION 

ThST FUk M = last COLU^i'm 

1 *in IF (M-M) I 5SM ^JO, I iS 
13S = 

GU TO nn 

TFST FOk L = StCOrjD LAST COLUMN 

l^iO IFCl-(N-D) igSMSo,ias 
liiS L = L*1 

GO TU SS 

ISO IF(If'JD-l) 160, iSSMbO 
ISS l<Nl() = 0 

GO T(t SO 

aV'^kAWt THPESmOLO aITFj final NOK'M 
160 IFCTHh^-ANk-^X) 16S, l6S#tiS 

SUkT hlGENVALDFS AND tIGFNVtCTOkS 


16S ID=-N 


PAGb 0003 


FIGEl 170 
F IGF 1 IPO 
F. IGF I 190 
FiGfc 1 poo 
FIGF. 1 PI 0 
F IGF 1 PPO 
E IGb 1 P30 
F IGF IPUO 
ElfiF 1 PSO 
E IGK1P60 
EIGF 1P70 
t 1(«F 1 kko 
E IGKlP'^o 
eige moo 

EK;F 1 SI 0 
EIGF 1 SPO 
EIGE 1330 
E IGE 1 Aun 
ElGEliSO 
F lliF 1 360 
EIGF 1 370 
EIGF 1 3H0 
EIGF. 1 390 
F ir;E M400 
E IGF 1 ai 0 
RIGFIUPO 
EIGE laxo 
EIGF 1 aoo 
EIGF 1 usn 
EIGF 1^60 
E IGE 1 070 
F I r; E W4 H 0 
e 1 GF 1 090 
F IGF 1 SOO 
E.I‘,F ISIO 
EIGF ISPO 
ElGh 1S30 
EIGF iS^iO 
F IGF ISSO 
Eir;E 1 ShO 
F 1GE1S70 
EIGE ISHO 
E1GFJS90 
EIGE1600 
EIGE 1610 
ETGE16P0 
EIGE 1630 
EIGElsao 
E1GE16S0 
E1GE1660 
E1GE1670 
EIGF 16H0 
EIGKl 690 
E1GE1700 
EIGF 1710 
EIGE 17P0 
EIGE1730 
EIGEi7an 


OI-V 


MEMHEfi Eir,EN 

00 1 = U‘M 

LL = i ♦( 1 ^n/^ 

= t J -rM 
O') 1>*S J=l,N 

v^v.rj^(j*j-J)/;> 

IF ( A(LL)-A(V^M) ) WO, l«S# 1 HS 
1 70 X = A(LL J 
A(LL 
Af''''5=x 

U (‘ V-1 ) WS, J 1 75 
WS 0 } 1 ^0 \ = 1 

1 L = 1 u< ♦ K 
I -Iwz J-J 4 «< 

x=M»nL^) 

■ 7 ( II W)=^( ) 

IPO y ( I ) =X 

IHb C )M I\Ut 

PE 

F.OO 


EIGEWSO 
EIGE 1 760 
ETGE W70 
tlGK 1 7P0 
tI(;tW90 
F IGF. WOO 
£Hit 1»510 
eu;f W<^0 
f:k;e w^o 

F. IGF wao 
F.lGtl HSO 
FlGKWhO 
t n, E W 7 n 
t IGhl.OPO 
fclGF. IFvOO 
F. IGF. WOO 

fk;f. WIO 

EJGMOi>n 
tIGE WAO 
EIGt jogo 
fc IGE WbO 


ooooor>oor:r>or!or>oononoooor^or»oonorir>onoonoooor>oooor>rjoonr>r)r)i 




C 

C 


NA.Vf DMF so 


SUHPOUTliJh 0'^‘FSO 
PUf<P(KSE 

FACTUW A t;iVfN SYV.SFTRIC POSITIVE OEFINITF MATRIX 


USA(,F. 

CALL 0 MF 50 C AfNrtPS, ItR) 


l-'tSCRIRnti^J 

A 


t PS 
lEW 


Of- PAkAf'FTEWS 

DiliJ^LF P^f-LISION UPPhW TRlArJOULAR PART OF UIVFN 

Sy'^v’KTRIL RUSUlVh DEFINITE :vl Hy N1 COE^’FICIENT 

MATRIX. 

i)N RFlJRN A CONTAI-jS THF RpSULTAMT UPPER 

tria,v«;ular matrix i<m oouhlf precision. 

Thk >jUMMeR ()P ^0.%S (COLn-'NSl IN :;iVEN matrix. 

SINOLF PrECISIOO I'jPtJT Cf'NSrANT WHICH IS USFO 

AS R&LATlvE tolerance FOR TFST ON LOSS OF 

blOr.JF ICAOCE. 

-feS'lLTINr; ERROR PARAvfTFR rODFO AS FOLLOWS 

ItPzO - iV') fhror 

ItRr-1 - Nl) RESULT RFCAOSF OF WWfiNO JNPlJT PARANE- 
TFR N OR HFCA1ISF sr)VE PAr>lC:A\0 IS NON- 
PDSiriVt (MATRIX ^ IS NOT POSITIVE 
OfFlom, possibly OIIF Tn LOSS OF SIONI- 
F I C Af^LE ) 

IE» = K - AARNlNf, which INDKATFS ( f)SS OF SK;N1FI- 
CA'iO.. THE RAOICArjO FH^viFO AT FACTORT/A- 
IION STEP K4.*i ^r.AS STILL POSITIVE hoT NO 
LONSFR OREATFR Ikan AhS C FPS* A ( h f 1 r k ♦ 1 ) ) . 


DMSD 

OMSO 

OMSO 

OMSO 

0 M 5 D 

OMsn 

DMSD 

DMSD 

DMSD 

DMSO 

OMSO 

()M$0 

DMSf) 

O'lSD 

DMSD 

DMSD 

OMSO 

DMSD 

DMSD 

DMSI> 

OMSO 

DMSD 

DMSD 

DMSD 

DMSD 

DffSI) 

DMSD 

DOSl) 

fJMSO 

DMSD 

DMSD 

DMSD 


P F M A H K S 

IHF UPPER IPIAMOOLAR PART QF OlVF.O MATRIX TS ASSUMED TO HF 
STORED COLIlMgwISF IN o(N + n/P SULCFSStVE STORAGE LOCATIONS, 

IN iHt same storage locatii)ns thf rfsoltin*; upper iwi angu- 
lar matrix is SKJRr.D CoLUMRi/TSF TOO. 

THF PRuCEDURF <,IVES RESULTS IF N 1$ GREATFR THAN 0 AND ALL 
CALCULATED RADICAN-:)5 ArE POSITIVE* 

TmF PRuDUCI of REUIRNEO DIAGONAL TERMS IS F.OUAL TO IhE 
SO'JARE-ROOT of the DFIERMIOANT of The GIVFN MATRIX; 

suHRCNiTiNFs and Function suhprograms rfooireo 

0 Li *\ E 
MF ThOO 

SOLUTlOiJ IS done using THE SO' I ARE -ROD T Mf rnOD OF CHOLESKY. 

TWt- r. TV^t^J MATO TV !C ^ t U ^ L c C k. T l: rv a o .1 .x . . rxc t <>x -r . x 1 , , r 


,DMSf) 

DMSD 

DMSD 

DMSD 

DMSD 

DMSD 


SUPRUUT IMF i)MF SIX A , N, EPS » J HR ) 

DIimENSION A(1) 


1 0 
?0 
30 
^JO 
SO 
80 


70 
80 
90 
100 
1 1 0 
1?0 
130 

IRO 
ISO 
180 
170 
180 
190 
POO 
?1 0 
PPO 
P^o 
Pu 0 
PSO 
880 
870 
?H 0 
890 


XJO 

Mo 

<80 

^SO 

380 

370 

38 0 
3 R 0 
iiOO 

ai 0 
p?o 
0 io 
aao 

^SO 

a80 

^70 

^80 

a9o 

soo 

SI 0 
S80 
S30 
Sao 

SSO 

S 80 

S70 

S 80 
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or» on ooo on on on no 


PARB 0006 


M£MRLft 

C 

C 


? 

3 

a 

s 

6 

7 

8 

9 


C 

C 


c 

c 


to 

u 


1? 


>JAVf: OMFSO 

[V.UJPLE PPeCISlON 0PIVrDSUvi,A 

IF ST ON .rPONf; INPUT PAPftMtftP N 
IF (N -1 ) iPr 1 r 1 

IMtlALUi 01 AUMNAL-LOOP 
KPIV20 

on 11 K = i 

•<PJ\/ = nPlVt‘^ 

T'vU = «PlV 
L t N 0 = K - 1 

CALCULATt TOLhPArJCF 
Ti.JL = A-‘S(tPS*SNt.;L( A(KPI V) ) ) 

STAkT FACTOkIZAl 1UN-L0(1P OVfcR K-TH POrt 
on 11 I=K,N 
OSlJ^'rO .00 
IF (I t g,)) 

STA»VT iN.JtP LOOP 
no 3 L = 1 f LF.ND 
LAN^-=KPI V-L 
LlNn=lNi>L 

DSU^' = nSU^^ + A (LANF)*fl CLINO) 
tfvO OF INf,EM LOOP 

•1PANSF0P*-1 hLFME.MT A(IND) 
nSU-*= A C INO) -USU^ 

IF ( I-K ) lOrS# 1 0 

Ttsi FUP NtGATIVfc PIVOT FLtM£NT ANM) FOR LOSS OF S I GN I F ] C A NC E 
IF (SNGL (lJ3U'i)-Tl)L) brS/9 
IF lO.SU *’) IPr 1^# 7 
IF ( ItP) 9 

IbP=K-l 

CUMPilTt PIVOT ELEMENT 
OPI VsDSOPT fOSUM) 

A(KPlV)=nPl V 
03! V=1 .00/UPlV 
GO T(‘ 11 

calculate terms in RO'a 
A( INiM=nSUM*OPl V 
lM)sl NOf] 

f NtJ OF Ul ARONAL-LUOP 

RFTUPNi 

It^=-1 

PETOh.'j 

FnO 


DMSD S90 
OMSD 60n 
OMSD 61 0 
oMsn 6 ?o 
OMSO 630 
OMSD 6^i0 
OMSD 6S0 
DMbO 660 
OMSD 670 
OMSD 6M0 
OMSD 690 
L)MSn 700 
OMSO 710 
OMSD 7P0 
DMSO 730 
OMSO 7ao 
OMSD 7S0 
l)MSO 760 
OMSO 770 
DMSn 7H0 
0'':.SO 790 
OMSD 800 
OMSD 810 
OMSU 820 
DMSO 830 
D^4SD 890 
DMSO 8S0 
OMSD 860 
DMSO 870 
OMSD 

DMSO 890 
DMSO 900 
DMSO 910 
DMSl) 9P0 
DMSO 930 
OMSD 990 
OMSD 9S0 
DMSO 960 
DMSO 970 
OMSD 9«0 
D'^'SD 990 
DMSDl OOn 
DMSOl 01 0 
OvsiM npo 
DMSOl 030 
DMSOl 090 
DMSOl OSO 
DMSDl 060 
OMS01070 
DMSOl 080 
DMSOl 090 
DMSDl lOQ 
DMSOl 110 
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PA(,h 0007 


vt‘-’RE W 
C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

f 

c 

c 

c 

c 

c 

c 

r 

c 

c 

c 

c 

c 

c 

c 

c 

c. 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

t: 

c. 

c 

r 

c 

c 


c 

c 

c 


NAMt DSTNV 


r l‘>jf DSINV 

p { 1 u tj n c p 

INvUt a GlVEfJ SY^METWIC POSITIVE np.FlNlTF MATWIX 
I'SAGt 

CALL OS I f'JV ( A , M , fc't^S f I LW ) 

Of-SCwlPTlO^^ C UPPKP TWUnGULAP ^ J c c t r t J 

SYwptPIC POSITIVE DEFlOITp N rtV N COEFFICIENT 

O'VrMU^^I a rO'^TAINS THF. WFStll. t ANT llFPt R 
tKlA.-KitiLAU MATRIX IN OOltHLE I S I ON , 

}, - iMt N-J-'HrR Of- RU'AS (COLl)MN.Sl IN r.IVIN 

FPS - SINCLF f’PtClSltlN INPUT CONSTANT *iHlCH IS U.StD 

AS PKUAIIXF TOLfcRANCt EOP TEST ON LOSS OF 

It.p - KtSlill 1 1 Ng'^'eRROP PARAMfcTPR COOFO AS FOLLOWS 


1 FR = 0 


Nci RFsI'LT UFCAI.ISE OE ^RONG INPUT PARAMfc- 
;j iiR HTCANSE S0*'E RADlCAM) IS NON- 
POGIIIVfc (VAIWIX A IS NOT fOSITIVe 

oFH-Niifc. pdssinly due to loss of SIGNI- 


POST I I Vfc 
OFF I-Nl It 
F U ANCt ) 
aAP ilNU 


ItRsK - aARnIOI. ,'<HICH IMOirATFS l.OSS OF SU'.NJFI- 
CAMCt. ImE RAOICAnO FORMFI) at FACIORI/A- 
TIIIN STEP K + 1 -.AS STILL POSITlVF ROT NO 
LONGtR GREATER Iman ABS ( ERi>* 1 • 

‘^^’"thf^upper triangular part of given matrix is *ssumeo to oe 

SlUPtO COLUVNrtlSF IN N.(N+1)/? SUCCESSIVE STORAGE LOCATIOWS, 
J,, Jet SAMt STORAGE inCATlOuS THE RESULTING UPPER IRI.ANGU- 
LAR 'lAIPlX IS STOffO CULU-'N‘”ISE TOIJ. 

IHE RWOCFUUPF GIVES RFSULIS IF N IS GREAIFP THAN 0 AND ALL 
CALCULATtU RAOICANOS ARE POSITIVE. 

SUOhOOTI Its AND FUNCTION SUHPHOGRAms RFQUIRfO 
O-'FSO 

^SOLOI lUN IS DONE USING F AC TOR 1 Z A I I ON BY SUBROUTINE OMFSI). 


S'JRRfUiT INt rSlNV ( A ,N,EP.S. IFR ) 


OIVENSIOM A(n 

DOUHLF PRECISION A,0IN,ft()RK 


FACIORl/E GlVt-N MATRIX hY '-’fANS OF SUliROUTlNF OMFSD 
A = IpAijSPUSE ( T ) • T 

CALL 0.••'FS0(A,..,EPS, IE») 


OSIN <3h0 
OSIN S70 
DSIN SSO 


^I-V 



ooo 



>JAM5 ^SI^V 
IFflh^O 9,1,1 


C 

C 


C 

C 


C 

C 


C 

C 


C 

C 


C 

C 

C 

c 

c 


c 

c 


L 

C 


c 

c 


Tvivt^r ijpppR r^^iAi\r,iiL AW vatwix t 
. A^yt T - JVt W$I MiNj-LnuP 
1 IwjVs'j* + l )/? 
ir.|f)=IPlV • ■ 

I*it UALWt iMVlr WSION-LOOP 
f)fi f> I = t , ^ 

Ol'jrl .(jO/A( IPTV) 

A ( 1 H I V ) = 0 I .‘j 

^ ft ' ; f z < • ^ j 

(-tw.)J S,S,? 

P .1 = 1 a) 

I ji riALi/b wo/^-i.nop 
00 u < = i,Ktr^u 
rw K z a , n 0 

L nny r 1 P I V 
LV£h=J . 

S1:kT LOOP 

00 5 L=LArjF,vi‘N 

LVb-5 = LVrW*l 

LH*jP = L»-cw-fL 

5 ^0W^ = A0^^KfA (LVFP)*A(LH0«) 
i) OF iNOtW LOOP 


A(J)=-l-.OPK*niN 
u j = j-'**i’g 

F. ‘‘O Or ^lJv\-LnOP 

S IPlv = TP;v-^:I;^l 
b rji)= 100-1 

F'jn OF iNvfcKSlOO-LOOP 


CALCOL ATK 
lO^FFSt (A) 

lOl T T AL 1/t 

00 H I r 1 , ij 

IPI V=IPI V^I 
.J= M*1 V 


lNV£PvSt(A) HY M£A\S OF INVr«SF(n 
= INVtPSF:(T) * T PANSPOSF. ( IMVF PSF ( n ) 
^ULT IPLICAT1UN-L0()P 


loi r I ALlZt PO>rj-LOOP 
OO ^ K = J , o 
\fjpK = o.oo 
LHllP=J 


STAkT IAimFM loop 
00 7 L = K,IJ 

LVFP = LMi)PfK_i 

AOKK = ..i.‘r^K ♦ A ( LHtiW ) * A (LVF9 ) 
7 L^U)»- =L-i(*.‘f -*-L 

E.W ) OF lOOtW LOOP 


PAoe: oooe 


DSIN S90 
OSIN 600 
OSIN 610 
OSIN 6P0 
OSIN 6iO 
0 SIN h a 0 
orufj 6S0 
OSIN 660 
OSIN h70 
OSIN 6 H n 
OSIN h90 
OMN 700 
OMN 710 
OSIN 7f^{) 
0-SIN 7.^0 
OSIN 7U0 
I5SI.N /SO 
OSIN 7hO 
OSIN 770 
OSIN 780 
OSIN 7«0 
OSIN 800 
OSI.M 810 
OSIN H^O 
OSIN 8 SO 
OSIN san 
OSIN 860 

osirg 8f7 0 
OSIN 870 
OSIN 860 
OSIN 890 
OSIN 900 
OSIN 910 
OSIN 980 
OSIN 9 So 
OSIN 9/1 f) 
OSIN 960 
OSJM 9/^0 
OSIN 970 
OSIN 980 
OSIN 990 
OSINJ 000 
OSIN 1010 
OS INI 080 
OS INI 0^0 
OSiNlOao 
OMNIOSO 
OS J N1 060 
ns INI 0 70 
OSTNl 080 
OSIN1090 
DSIN 1 1 00 
OSINl 1 10 
OSIM 1?0 
OSINl 1 ^0 
OSINl 1 ao 
OSINl ISO 
OSINl 160 
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PAUE 0009 



lEMHfP NflMK nSI<'JV 
A ( J ) = 0 w K 
b J=J>K 

C UF klj/.- AND MULTIPlICATIUN-LOOF^ 

c 

9 

t M •' J 


DSINl 170 

oslNlibn 

OSlfvM 190 
I>SUj1?00 
05IN1 ;?1 0 
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